home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / atan2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  5.7 KB  |  281 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *   Replacement atan2() function for PML library. Original function    *
  4.  *   from Fred Fish does not conform to Standard C.  It has arguments   *
  5.  *   in a reverse order.  It also does not complain if both arguments   *
  6.  *   are 0.0                                                            *
  7.  *                                    *
  8.  *   Michal Jaegermann - December 1990
  9.  *                                    *
  10.  ************************************************************************
  11.  */
  12. /************************************************************************
  13.  *                                    *
  14.  *                N O T I C E                *
  15.  *                                    *
  16.  *            Copyright Abandoned, 1987, Fred Fish        *
  17.  *                                    *
  18.  *    This previously copyrighted work has been placed into the    *
  19.  *    public domain by the author (Fred Fish) and may be freely used    *
  20.  *    for any purpose, private or commercial.  I would appreciate    *
  21.  *    it, as a courtesy, if this notice is left in all copies and    *
  22.  *    derivative works.  Thank you, and enjoy...            *
  23.  *                                    *
  24.  *    The author makes no warranty of any kind with respect to this    *
  25.  *    product and explicitly disclaims any implied warranties of    *
  26.  *    merchantability or fitness for any particular purpose.        *
  27.  *                                    *
  28.  ************************************************************************
  29.  */
  30.  
  31.  
  32. /*
  33.  *  FUNCTION
  34.  *
  35.  *    atan2   double precision arc tangent of two arguments
  36.  *
  37.  *  KEY WORDS
  38.  *
  39.  *    atan2
  40.  *    machine independent routines
  41.  *    trigonometric functions
  42.  *    math libraries
  43.  *
  44.  *  DESCRIPTION
  45.  *
  46.  *    Returns double precision arc tangent of two
  47.  *    double precision floating point arguments ( atan(Y/X) ).
  48.  *
  49.  *  USAGE
  50.  *
  51.  *    double atan2(y, x)
  52.  *    double y;
  53.  *    double x;
  54.  *
  55.  *  (This order of arguments really makes sense if you will notice
  56.  *   that this function is mainly used during a conversion from
  57.  *   rectangular to polar coordinates)
  58.  *
  59.  *  REFERENCES
  60.  *
  61.  *    Fortran 77 user's guide, Digital Equipment Corp. pp B-4.
  62.  *
  63.  *  RESTRICTIONS
  64.  *
  65.  *    For precision information refer to documentation of the
  66.  *    other floating point library routines called.
  67.  *    
  68.  *  PROGRAMMER
  69.  *
  70.  *    Fred Fish
  71.  *    Modified by Michal Jaegermann
  72.  *
  73.  *  INTERNALS
  74.  *
  75.  *    Computes atan2(y, x) from:
  76.  *
  77.  *        1.    If x = 0 and y != 0 then
  78.  *            atan2(y, x) = PI/2 * (sign(y))
  79.  *                      otherwise error signaled and
  80.  *                      atan2(y ,x) is 0.
  81.  *
  82.  *        2.    If x > 0 then
  83.  *            atan2(y, x) = atan(y/x)
  84.  *
  85.  *        3.    If x < 0 then atan2(y, x) =
  86.  *            PI*(sign(y)) + atan(y/x)
  87.  *
  88.  */
  89.  
  90. #if !defined (__M68881__) && !defined (sfp004)
  91.  
  92. #include <stdio.h>
  93. #include <math.h>
  94. #include "pml.h"
  95.  
  96. static char funcname[] = "atan2";
  97.  
  98.     struct exception xcpt;
  99.  
  100. double atan2 (y, x)
  101. double y;
  102. double x;
  103. {
  104.     double result;
  105.  
  106.     if (x == 0.0) {
  107.     if (y == 0.0) {
  108.         result = HUGE_VAL;    /* mjr: for now, should be NAN    */
  109.         xcpt.type = DOMAIN;
  110.         xcpt.name = funcname;
  111.         xcpt.arg1 = x;
  112.         xcpt.retval = result;
  113.         if (!matherr(&xcpt)) {
  114.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  115.         errno = EDOM;
  116.         }
  117.     } else {
  118.         result = copysign (HALFPI, y);
  119.     }
  120.     } else {
  121.     result = atan (y/x);
  122.     if (x < 0.0) {
  123.         result += copysign (PI, y);
  124.     }
  125.     }
  126.     return (result);
  127. }
  128.  
  129. #endif !defined (__M68881__) && !defined (sfp004)
  130. #ifdef    __M68881__
  131. __asm("
  132. .text
  133. .even
  134. _funcname:
  135.     .ascii    \"atan2\\0\"
  136.     .even
  137.  
  138. .globl    _atan2
  139. _atan2:
  140. | denormalized numbers are treated as 0
  141.     tstl    sp@(12)
  142.     beq    5f        | x == 0!
  143.     blt    1f        | x < 0!
  144.                 | x > 0: return atan(y/x)
  145.  
  146.     fmoved    sp@(4),fp0    | get y
  147.     fdivd    sp@(12),fp0    | y/x    
  148.     fatanx    fp0,fp0        | atan(y/x)
  149.     bra 3f            | return
  150. 1:                | x < 0
  151.  
  152.     fmovecr    #0,fp1        | get pi
  153.     fmoved    sp@(4),fp0    | get y
  154.     fdivd    sp@(12),fp0    | y/x
  155.     fatanx    fp0,fp0        | atan(y/x)
  156.     btst    #31,sp@(4)    | sign(y)
  157.     beq    2f        | positive!
  158.  
  159.     fnegx    fp1,fp1        | transfer sign
  160. 2:    faddx    fp1,fp0        | sign(y)*pi + atan(y/x)
  161. |    bra 3f            | return
  162. 3:
  163.     fmoved    fp0,sp@-    | return result
  164.     moveml    sp@+,d0/d1
  165. 4:    
  166.     rts            | sigh.
  167. 5:                | x == 0
  168.     movel    #1073291771,d0    | pi/2
  169.     movel    #1413754136,d1    |
  170.  
  171.     tstl    sp@(4)        | 
  172.     beq    6f        | NaN
  173.     bge    4b        | exit
  174.     bset    #31,d0        | x < 0 : return -pi/2
  175.     bra    4b
  176. 6:    movel    #-1,d0        | NaN
  177.     movel    #-1,d1        |
  178.     bra    4b
  179. ");    /* end asm    */
  180. #endif    __M68881__
  181.  
  182. #ifdef    sfp004
  183. __asm("
  184.  
  185. comm =     -6
  186. resp =    -16
  187. zahl =      0
  188.  
  189. .even
  190. .text
  191. .even
  192. _funcname:
  193.     .ascii    \"atan2\\0\"
  194.     .even
  195. .text
  196. .even
  197. .globl    _atan2
  198. _atan2:
  199. | denormalized numbers are treated as 0
  200.     lea    0xfffa50,a0
  201.     moveml    a7@(12),d0-d1    |  x
  202.     tstl    d0
  203.     beq    5f        | x == 0!
  204.     blt    1f        | x < 0!
  205.                 | x > 0: return atan(y/x)
  206.  
  207. |    fmoved    sp@(4),fp0    | get y
  208.     movew    #0x5400,a0@(comm)
  209.     .long    0x0c688900, 0xfff067f8
  210.     movel    sp@(4),a0@
  211.     movel    sp@(8),a0@
  212.  
  213. |    fdivd    sp@(12),fp0    | y/x
  214.     movew    #0x5420,a0@(comm)
  215.     .long    0x0c688900, 0xfff067f8
  216.     movel    d0,a0@
  217.     movel    d1,a0@
  218.  
  219. |    fatanx    fp0,fp0        | atan(y/x)
  220.     movew    #0x000a,a0@(comm)
  221.     .word    0x4a68,0xfff0,0x6bfa
  222.  
  223.     bra 3f            | return
  224. 1:                | x < 0
  225.  
  226. |    fmovecr    #0,fp1        | get pi
  227.     movew    #0x5c80,a0@(comm)
  228.     .long    0x0c688900, 0xfff067f8
  229.  
  230. |    fmoved    sp@(4),fp0    | get y
  231.     movew    #0x5400,a0@(comm)
  232.     .long    0x0c688900, 0xfff067f8
  233.     movel    sp@(4),a0@
  234.     movel    sp@(8),a0@
  235.  
  236. |    fdivd    sp@(12),fp0    | y/x
  237.     movew    #0x5420,a0@(comm)
  238.     .long    0x0c688900, 0xfff067f8
  239.     movel    d0,a0@
  240.     movel    d1,a0@
  241.  
  242. |    fatanx    fp0,fp0        | atan(y/x)
  243.     movew    #0x000a,a0@(comm)
  244.     .word    0x4a68,0xfff0,0x6bfa
  245.  
  246.     btst    #31,sp@(4)    | sign(y)
  247.     beq    2f        | positive!
  248.  
  249. |    fnegx    fp1,fp1        | transfer sign
  250.     movew    #0x049a,a0@(comm)
  251.     .word    0x4a68,0xfff0,0x6bfa
  252.  
  253. 2:|    faddx    fp1,fp0        | sign(y)*pi + atan(y/x)
  254.     movew    #0x0422,a0@(comm)
  255.     .word    0x4a68,0xfff0,0x6bfa
  256.  
  257. |    bra 3f            | return
  258. 3:
  259. |    fmoved    fp0,d0-d1    | return result
  260.     movew    #0x7400,a0@(comm)
  261.     .long    0x0c688900, 0xfff067f8
  262.     movel    a0@,d0
  263.     movel    a0@,d1
  264.  
  265. 4:    
  266.     rts            | sigh.
  267. 5:                | x == 0
  268.     movel    #1073291771,d0    | pi/2
  269.     movel    #1413754136,d1    |
  270.  
  271.     tstl    sp@(4)        | 
  272.     beq    6f        | NaN
  273.     bge    4b        | exit
  274.     bset    #31,d0        | x < 0 : return -pi/2
  275.     bra    4b
  276. 6:    movel    #-1,d0        | NaN
  277.     movel    #-1,d1        |
  278.     bra    4b
  279. ");    /* end asm    */
  280. #endif    sfp004
  281.